library(readxl)
library(tidyverse)
library(plyr) # mapvalues()
library(descr) # freq()
library(DescTools) # statistic tests
library(modelsummary)
library(marginaleffects)
library(gridExtra) # combine graphs
library(lemon) # put legend in combined graphs
library(Hmisc) # correlation matrix & weighted sd
library(DescTools) # pseudo r-squared

PRRdataframe1 <- read_excel("Dataset_rfj.xlsx")
PRRdataframe2 <- na.omit(PRRdataframe1)

freq(PRRdataframe2$Party)

8807/2
4403.5/887 # 4.96
4403.5/7920 # 0.56

PRRdataframe2$df_weights <- mapvalues(PRRdataframe2$Party,
                                      from = c(0, 1),
                                      to = c(4.96, 0.56))

### PCAs on activism ###

# 1st run

PCA_dataframe <- PRRdataframe2[, c(13:17)]
pca_result <- prcomp(PCA_dataframe, scale = TRUE)
pca_result$rotation <- -1*pca_result$rotation
pca_result$x <- -1*pca_result$x

summary(pca_result)

# Scree plot
var_explained = pca_result$sdev^2 / sum(pca_result$sdev^2)

qplot(c(1:5), var_explained) + 
  geom_line() + 
  xlab("Principal Component") + 
  ylab("Variance Explained") +
  ggtitle("Scree Plot") +
  ylim(0, 1) +
  theme_bw()

# 2nd run

PCA_dataframe <- PRRdataframe2[, c(13, 15:17)]
pca_result <- prcomp(PCA_dataframe, scale = TRUE)
pca_result$rotation <- -1*pca_result$rotation
pca_result$x <- -1*pca_result$x

summary(pca_result)

# Scree plot
var_explained = pca_result$sdev^2 / sum(pca_result$sdev^2)

qplot(c(1:4), var_explained) + 
  geom_line() + 
  xlab("Principal Component") + 
  ylab("Variance Explained") +
  ggtitle("Scree plot") +
  ylim(0, 1) +
  theme_bw()

PCA_activism <- (pca_result$x[, 1]) # Extract scores of the first principal component (PC1)
PRRdataframe2$PCA_activism <- PCA_activism # Add the new variable to your original data frame

#=========================================#

# Table

CrossTable(PRRdataframe2$Purposive, PRRdataframe2$Gender, prop.c=TRUE,
           prop.t=FALSE, prop.r=FALSE, prop.chisq = FALSE,
           total.r=FALSE, total.c=TRUE, chisq=TRUE)

CrossTable(PRRdataframe2$Solidary, PRRdataframe2$Gender, prop.c=TRUE,
           prop.t=FALSE, prop.r=FALSE, prop.chisq = FALSE,
           total.r=FALSE, total.c=TRUE, chisq=TRUE)

CrossTable(PRRdataframe2$Material, PRRdataframe2$Gender, prop.c=TRUE,
           prop.t=FALSE, prop.r=FALSE, prop.chisq = FALSE,
           total.r=FALSE, total.c=TRUE, chisq=TRUE)

CrossTable(PRRdataframe2$Expressive, PRRdataframe2$Gender, prop.c=TRUE,
           prop.t=FALSE, prop.r=FALSE, prop.chisq = FALSE,
           total.r=FALSE, total.c=TRUE, chisq=TRUE)

CrossTable(PRRdataframe2$Altruistic, PRRdataframe2$Gender, prop.c=TRUE,
           prop.t=FALSE, prop.r=FALSE, prop.chisq = FALSE,
           total.r=FALSE, total.c=TRUE, chisq=TRUE)

#=========================================#

# Model

coefficientsnames <- c("Party" = "Party (SD)",
                       "yearsmember" = "Years of membership",
                       "Labourforce" = "Active in the workforce",
                       "Tertiaryedu" = "Tertiary education",
                       "agesquared" = "Age squared",
                       "age" = "Age",
                       "Altruistic" = "Group incentives",
                       "Expressive" = "Expressive incentives",
                       "Material" = "Material incentives",
                       "Solidary" = "Solidary incentives",
                       "Purposive" = "Purposive incentives",
                       "Gender" = "Gender (Woman)")

PCA_fullmodel <- lm(PCA_activism ~ Gender + Purposive + Solidary + Material + Expressive + Altruistic +
             age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
           data = PRRdataframe2, weights = df_weights)
modelsummary(PCA_fullmodel, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
modelsummary(PCA_fullmodel, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "PCA full model.docx")

# Material incentives

activism_0_material <- lm(PCA_activism ~ Gender + Purposive + Solidary + Expressive + Altruistic +
                            age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                         data = PRRdataframe2, weights = df_weights)
modelsummary(activism_0_material, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

activism_M_material <- glm(Material ~ Gender + Purposive + Solidary + Expressive + Altruistic +
                             age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                          data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(activism_M_material, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(activism_M_material)

# Group incentives

activism_0_group <- lm(PCA_activism ~ Gender + Purposive + Solidary + Material + Expressive +
                            age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                          data = PRRdataframe2, weights = df_weights)
modelsummary(activism_0_group, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

activism_M_group <- glm(Altruistic ~ Gender + Purposive + Solidary + Material + Expressive +
                             age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                           data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(activism_M_group, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(activism_M_material)


